On a vu en M1 ENSP :
Source : allison
Horst @allison_horst
Objectifs en M2 ENSP :
Source : allison
Horst @allison_horst
Toujours les mêmes bonnes habitudes de travail :
-> Exprimez vos difficultés et blocages : ce qui permet de reprendre ensemble des aspects vus trop rapidement et d’éviter que les difficultés s’accumulent jusqu’à devenir insurmontables
-> Jouez le jeu des exercices : c’est un moyen (le seul) de se familiariser avec l’écriture de R : Les réponses sont sur les supports de cours mais faites l’effort de chercher par vous-même.
Est-ce que ce programme vous convient et est complémentaire d’autres cours ? Pas de doublons ? D’autres attentes ?
Les différents jeux de données utilisés en cours seront présentés au fur et à mesure.
Il s’agira d’un dossier en groupe de 2 à 4 étudiant·es au choix à partir de données fournies et à rechercher autour d’une même thématique. Le dossier rendu sera :
-> Vous m’enverrez le document produit et le script .Rmd au plus tard le 03/01/23 à mon adresse mail cecile.rodrigues@univ-lille.fr (vous aurez un accusé de réception sous 24h).
Objectif : Revoir tout le cours de M1 en 3 heures
R est un dérivé gratuit (licence GNU GPL) du langage de statistique S. C’est un projet collaboratif et donc ouvert à tous. Il est entretenu par le R Development Core Team (CRAN). C’est un langage permettant de produire des statistiques descriptives comme des statistiques avancées. Ses fonctionnalités sont étendues ; il offre un large choix d’options graphiques. Il permet de développer ses propres fonctions et d’utiliser des fonctions développées par d’autres chercheurs, ingénieurs ou développeurs.
En bref :
… Quelques inconvénients :
La console disponible avec R (RGui) est rudimentaire. Il est conseillé d’utiliser un logiciel offrant un espace de travail plus agréable. RStudio est un environnement de développement pour l’écriture et la visualisation de scripts R. Le logiciel est produit par une société (qui va prochainement s’appeler posit) à l’origine de nombreuses ressources et projets autour de R ({Rmarkdown} pour les rapports, {shiny} pour les applications, {tidyverse}, RstudioCloud, etc.). Il est également disponible gratuitement. Il offre un éditeur de texte (coloration, auto complétion, indentation) et un affichage simultané des scripts, sorties, fichiers, graphiques, aides. Il est quasiment systématiquement utilisé pour écrire en R.
- (1) Télécharger le langage R :
sudo apt-get install r-base r-base-dev- (2) Télécharger et installer l’IDE RStudio :
Attention :
Packages : Extensions ajoutées à votre session de travail, une fois installées et chargées. Les packages sont des bibliothèques de fonctions supplémentaires développées par d’autres personnes.
Les packages “officiels” sont validés par la communauté de R (CRAN). Ils sont tous accompagnés d’une documentation (document pdf issu du site cran.r-project.org).
Pour utiliser un package validé par le CRAN :
install.packages("tidyverse", dep=T)Cette étape d’installation des packages est délicate : R évoluant rapidement, la compatibilité entre les packages et entre les packages et la version utilisée de R n’est toujours maintenue.
library(tidyverse)
library(questionr)/!\ Il est préférable d’écrire ces lignes afin que les packages soient automatiquement définis et que les fonctions associées fonctionnent.
On recommande aussi de mettre ces lignes en début de script.
library() est équivalent à require()
Enregistrer votre script : File > Save as…
Quand on entame un projet, les lignes de commande, scripts, tableaux créés et dossiers peuvent s’accumuler. Il est important d’organiser son travail dès le début.
Gérer ses projets. Un projet est une session qu’on associe à un dossier :
Vous retrouvez votre session : script, objets, recodages et transformations déjà effectués
Renvoie directement à cet emplacement
Evite de mélanger les données et scripts portant sur différents sujets
Permet de partager l’ensemble du dossier à d’autres personnes
En haut à droite, une icône en forme de cube bleu permet de créer un “nouveau projet”. C’est un assistant de création d’un dossier dédié à un projet spécifique dans votre explorateur (windows, mac ou linux). Il sera associé à une session de Rstudio, qui sera indépendante des documents et objets utilisés pour un autre projet.
ou bien File > New Project > Existing directory –> Browse > Sélectionnez votre dossier créé pour le TD > Create project
Un fichier Nom_du_dossier.Rproj apparaît dans le dossier, c’est le raccourci qui mènera à la session R. Un fichier .RData sera également créé et contiendra les tableaux présents dans l’environnement.
Commenter ses scripts
Les commentaires sont précédés de # et ne sont pas exécutés.
Pour passer un ensemble sélectionné en commentaire (et inversement) : ctrl+shift+c.
Aérer le texte et mettre des titres facilitent aussi grandement la lecture de scripts : sur une ligne, # Titre ----
Dernier rappel
Il existe toujours plusieurs moyens d’aboutir à un résultat ; L’essentiel est de comprendre ce que vous recopiez / adaptez / écrivez et de favoriser les solutions les plus simples (pour faciliter la relecture, minimiser la mémoire sollicitée et les risques d’erreurs)
Veillez à :
Organiser vos différents scripts à l’aide des sections et commentaires : un·e futur·e vous ou un·e collègue doit pouvoir se replonger facilement dans le code
Ne laisser que le code utile ne comportant pas d’erreur ou d’instructions longues à exécuter et/ou inutile à exécuter (comme un install.packages). Alternez entre les zones console (retours et tests) et les script (code enregistré). Dans l’idéal, vous devriez pouvoir exécuter l’ensemble du script sans erreur (le bouton “source” exécute l’ensemble du script et s’arrête à la première erreur).
-> Si vous ne l’avez pas fait, créez un dossier consacré au TD de R avancé (M2), appelé par exemple “TD_R_M2”.
-> Téléchargez le support de cours pour le mettre dans ce dossier.
-> Créez un projet .Rproj associé au dossier du TD et ouvrez un script .R
Fichier > Nouveau projet > Existing directory > Browse > Choisir le dossier créé pour le TD
Avant toutes choses, commencez à lire la documentation associée aux premières données Coco1 qu’on utilisera dans les parties 1 à 3 du TD :
Documentation_coco1.pdf : pp. 4-8 ; La suite est le dictionnaire des codes
Questionnaire_coco1.pdf : L’ensemble ; En notant les codes des questions qui vous intéressent
!! Attention : Distinguez bien la Documentation qui comporte le dictionnaire des codes du Questionaire qui ne contient que les questions de l’enquête. Les données reçues contiennent en plus des réponses au questionnaire une série de variables socio-démographiques issues de l’enquête annuelle réalisée par l’équipe ELIPSS. Vous avez toujours le fichier Documentation_coco1.pdf ouvert pendant que vous travaillez sur les données Vous trouverez le questionnaire de l’enquête annuelle dans l’archive cdsp_ea2019, fichier Questionnaire_EA2019.pdf.
coco
de la vague 1Suivant les types des données à importer (voir leur extension), différentes fonctions d’importation seront utilisées :
readRDS("chemin/nomFichier.rds") ;read.csv2("chemin/nomFichier.csv", stringsAsFactors=F)
;read.table("chemin/nomFichier.txt", sep="\t")Un package utile facilitant les importations et exportations (entre
la session R et le dossier associé au projet) : rio : voir
la
vignette
de présentation
… /!\ : Beaucoup de packages facilitent les actions : risque de s’y perdre et risque d’être facilement bloqué·e en cas de bug si vous ne connaissez pas un minimum le R-base
unzip("data/cdsp_coco1_2020/cdsp_coco1_2020.zip", exdir="data/")
coco <- read.csv("data/cdsp_coco1_2020/donnees_coco1/csv_coco1/fr.cdsp.ddi.elipss.202004coco1.csv",
stringsAsFactors = F)- Repérez les noms des variables contenant le genre, l’âge et le niveau de diplôme des enquêté·es.
Le genre se trouve dans 2 variables (eayy_A1 et coco1_cal_SEXE), l’âge dans plusieurs (eayy_A2A_rec pour l’âge quinquennal et coco1_cal_AGE1 pour l’âge décennal) et le niveau de diplôme est dans la variable eayy_B18_rec.
La différence est que les variables préfixées eayy sont les variables issues de l’enquête annuelle (à utiliser) et les variables préfixées cal sont les variables de calages utilisées pour produire les coefficients de redressement et qu’on utilisera pas. On utilisera plutôt les variables eayy.
Effectifs avec table :
table(df$Variable)table(df$1èreVariable, df$2èmeVariable)Effectifs et poucentages sur une variable avec
freq : freq(df$variable)
Pourcentages avec prop.table(),
rprop() et cprop() s’appliquant sur une table
d’effectifs :
rprop(table(df$1èreVariable, df$2èmeVariable))
Indices de répartition : min,
max, sd, mean,
median et summary.
Exporter un résultat :
write.csv2(resultat, "resultats/tc_prepa_genre.csv") ;
s’assurer qu’il existe un dossier “resultat” puis mise en forme sous
Excel ou Calc
table$sexe[table$sexe=="1"] <- "Homme"# 2. Recodage :
table$sexe <- fct_recode(
factor(table$sexe),
"Homme"="1")
# ou fct_collapse()table$recodage[table$recodage=="Non renseigné"]<-NA
table$recodage[is.na(table$recodage)]<-"Non renseigné"
# Ou
table$recodage <- fct_recode(factor(table$variable_a_recoder),
"Valeur 1"="1","Valeur 2"="2",
NULL="Non renseigné")
table$recodage <- fct_explicit_na(table$recodage,
na_level = "Non renseigné")# 1. Chargement des packages utilisés :
library(questionr)
library(tidyverse)
# 2. Recodage :
coco$sexe <- fct_recode(
factor(coco$eayy_A1),
"Homme"="1", "Femme"="2")
# 3. Tableau de fréquences :
freq(coco$sexe,
total= T, valid=F)2 fonctions principales permettent de recoder des variables sur conditions :
ifelse(condition, si condition remplie, si condition non remplie)case_when( conditions 1 ~ effet de la condition 1 ,
condition 2 ~ effet de la condition 2,
etc,
TRUE ~ effet si aucune condition n'est remplie)# 1. Recodages :
coco$coco1_q19 <- fct_recode(factor(coco$coco1_q19),
"En permanence"="1",
"Souvent"="2",
"Quelques fois"="3",
"Rarement"="4",
"Jamais"="5",
NULL = "9999")
# 2. Tableau :
tc_genre_tristesse <- rprop(table(sexe=coco$sexe,
tristesse=coco$coco1_q19))
tc_genre_tristesse## tristesse
## sexe En permanence Souvent Quelques fois Rarement Jamais Total
## Homme 0.8 3.8 18.0 36.5 40.9 100.0
## Femme 1.1 9.5 32.8 29.6 27.0 100.0
## Ensemble 0.9 6.8 25.9 32.9 33.5 100.0
Source : allison
Horst @allison_horst
GGplot2 est un package dédié aux graphiques. Sa syntaxe
est cohérente avec celle de Tidyverse. L’intérêt de
Ggplot2 est de produire rapidement et efficacement des
graphiques complexes par ajout de couches successives. De plus, il est
possible d’enchaîner la mise en forme des tableaux avec la réalisation
des graphiques.
La syntaxe de base de ggplot est :
ggplot (table_de_données) + aes (x =
Variable_abscisses, y = Variable_ordonnées ) + geom_X
(elements de mise en forme) + labs (titres)
où ggplot() définit la source des données,
geom_X() le type de représentation, aes() les
données à représenter et labs() les différents titres et
textes.
A laquelle on peut ajouter des précisions du type :
theme() : Mise en forme des zones de texte et de
dessin
scale_y_continuous() : Options sur les
échelles
scale_fill_manual() : Options sur les
couleurs
# 4. Graphique :
tc_genre_tristesse %>%
data.frame() %>%
filter(tristesse !="Total") %>%
ggplot()+
aes(fill=sexe, x=tristesse, y= Freq) %>%
geom_bar(col="black",stat="identity", position="dodge")+
labs(title="Tristesse selon le genre pendant le 1er confinement",
x="S'est senti·e triste",y="%",fill="",
caption = "Source : Coco1 - ELIPSS, CDSP, avril 2020")+
viridis::scale_fill_viridis(discrete=T)+
theme_bw()On utilise les questions Q12 et Q13. Il y a peu de variables réellement numériques.
class(coco$coco1_q12)## [1] "integer"
coco$coco1_q12[coco$coco1_q12 %in% c(6666,9999)] <- NA
coco$coco1_q13[coco$coco1_q13 %in% c(6666,9999)] <- NA
coco$gueries <- round(coco$coco1_q13/coco$coco1_q12*100,1)
summary(coco$gueries)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 33.30 100.00 65.52 100.00 100.00 648
En fonction de l’âge ? En fonction de la PCS ?
%>%Les principales fonctions permettant de réarranger les tables via tidyverse sont :
select(table, colonnes à conserver) : Sélection de
colonnes ; Voir starts_with(), contains(),
ends_with()filter(table, condition du filtre) : Sélection de
lignesarrange(table, colonnes sur laquelle faire le tri, desc(2ème colonne sur laquelle faire un tri décroissant))
: Tri de la tableslice(table, position des lignes à conserver) :
Sélection de lignesmutate(table, nouvelleVariable = ) : Création de
variablesrename("nouveau nom" = "ancien nom") : Renommer une
variable%>% (pipe, ctrl+shift+m) permet
d’enchaîner des fonctions : f(a) s’écrit
a %>% f().Il peut être intéressant de regrouper plusieurs observations selon un
ou plusieurs critères afin de produire des calculs. La fonction
group_by permet de produire ces regroupements et est
associée à une fonction de calcul dans summarise.
Exemple : On souhaite compter le nombre d’individus par classe d’âge ainsi que le revenu moyen par classe d’âge d’une table donnée :
table %>%
group_by(classe_dage) %>%
summarise(Nb_par_classe_dage = n(),
Revenu_moyen = mean(revenu, na.rm=T))# 1. Recodage :
coco$age_quinq <- fct_recode(factor(coco$eayy_A2A_rec),
"<24 ans"="4", "25-29 ans"="5",
"30-34 ans"="6", "35-39 ans"="7",
"40-44 ans"="8", "45-49 ans"="9",
"50-54 ans"="10", "55-59 ans"="11",
"60-64 ans"="12","65-69 ans"="13",
"70 ans et+"="14")
# 2. Avec group_by :
coco %>%
group_by(age_quinq) %>%
summarise(mediane = median(gueries, na.rm=T),
moyenne = mean(gueries, na.rm=T))coco$pcs <- fct_collapse(factor(coco$eayy_PCS6CJT),
"Agriculteurs"="1",
"Artisans, comm., chef d'E"="2",
"Cadres et pr.intell.sup."="3",
"Prof. intermédiaires"="5",
"Employés"="4","Ouvriers"="6",
NULL = c("666","999"))
coco %>%
group_by(pcs) %>%
summarise(mediane = median(gueries, na.rm=T),
moyenne = mean(gueries, na.rm=T))names(table) permet d’agir sur les noms des
variables
gsub(pattern = "on", replacement = "a", x = "bonbon")
modifie tous les “on” en “a” de “bonbon”.
names(coco) <- gsub("coco1_","",names(coco))On passera par la création de variables intermédiaires :
sitpro.avt : Situation d’emploi avant le confinement
; A partir de la variable q32
sitpro.conf : Situation d’emploi pendant le
confinement ; A partir de la variable q37
tt.conf : Et on intègre le télétravail avec
q38
# freq(coco$q32)
coco$sitpro.avt <- fct_collapse(factor(coco$q32),
"En emploi, contrat" = c("1","2"),
"Sans emploi" = c("3","4"),
"Retraite" = "6",
"Autre inactif·ve" = c("5","7","8"), NULL="9999")
freq(coco$sitpro.avt, sort="dec", valid=F)# freq(coco$q37)
coco$sitpro.conf <- fct_collapse(factor(coco$q37),
"En congé"=c("1","2","3"),
"Au chômage"=c("4","5"),
"En emploi"=c("6"),
NULL=c("7","6666","9999"))
freq(coco$sitpro.conf)# freq(coco$q38)
coco$tt.conf <- fct_recode(factor(coco$q38),
"Lieu habituel"="1",
"Lieu habituel/TT"="2",
"Télétravail"="3", NULL = "6666")
freq(coco$tt.conf)# table(coco$sitpro.conf, coco$tt.conf, useNA = "ifany")
coco$sitpro.tt.conf <- ifelse(is.na(coco$tt.conf) ==F,
paste0("En emploi - ", coco$tt.conf),
as.character(coco$sitpro.conf))
freq(coco$sitpro.tt.conf)De fait, les questions sur leschangements de situation d’emploi ne concernent que les personnes précédemment en emploi. Pour analyser ces questions, on passe par la création d’une table conservant la sous-population des actif·ves avant la pandémie.
–
coco.activite <- coco %>% filter(sitpro.avt == "En emploi, contrat")
# freq(coco.activite$sitpro.tt.conf, sort="dec",valid=F)On peut maintenant envisager les modifications des conditions de travail selon l’âge, la PCS, etc.
rprop(table(coco.activite$pcs, coco.activite$sitpro.tt.conf,
useNA="ifany"))##
## Au chômage En congé En emploi
## Agriculteurs 25.0 12.5 0.0
## Artisans, comm., chef d'E 21.1 15.8 0.0
## Cadres et pr.intell.sup. 17.7 6.3 2.5
## Employés 18.0 15.3 7.2
## Prof. intermédiaires 24.8 20.0 3.8
## Ouvriers 29.0 14.5 4.3
## <NA> 25.6 10.6 6.0
## Ensemble 23.2 13.2 4.9
##
## En emploi - Lieu habituel
## Agriculteurs 50.0
## Artisans, comm., chef d'E 31.6
## Cadres et pr.intell.sup. 20.3
## Employés 19.8
## Prof. intermédiaires 29.5
## Ouvriers 30.4
## <NA> 21.6
## Ensemble 24.2
##
## En emploi - Lieu habituel/TT
## Agriculteurs 0.0
## Artisans, comm., chef d'E 10.5
## Cadres et pr.intell.sup. 6.3
## Employés 5.4
## Prof. intermédiaires 1.9
## Ouvriers 2.9
## <NA> 5.5
## Ensemble 4.7
##
## En emploi - Télétravail <NA> Total
## Agriculteurs 0.0 12.5 100.0
## Artisans, comm., chef d'E 21.1 0.0 100.0
## Cadres et pr.intell.sup. 45.6 1.3 100.0
## Employés 32.4 1.8 100.0
## Prof. intermédiaires 17.1 2.9 100.0
## Ouvriers 17.4 1.4 100.0
## <NA> 29.1 1.5 100.0
## Ensemble 27.8 1.9 100.0
!! Les effectifs commencent à devenir faibles (en lien avec la structure de la population du panel ELIPSS).
freq(coco.activite$pcs)Des catégories devraient être regroupées…
Pour en terminer avec les exercices de rappel du M1, on veut représenter graphiquement les risques professionnels encourus pendant la période de confinement pour les personnes déjà en activité en fonction du genre.
Les variables vont les q41_1 à q41_6 et viennent d’une question à choix multiples :
Si Q38 = {1, 2}
Q41 - Au cours des deux dernières semaines, vos conditions de travail vous ont-elles exposé à : (Plusieurs réponses possibles)
.. plutôt dommage qu’on ne pose pas la question aux personnes en télétravail.
On commence par les recodages :
# Recodage du risque pour la santé :
coco.activite$risque.sante <- fct_recode(factor(coco.activite$q41_1),
"Oui"="1","Non"="2",NULL = "6666")
freq(coco.activite$risque.sante)# Recodage du risque de conflit avec l'employeur :
coco.activite$risque.conflit.empl <- fct_recode(factor(coco.activite$q41_2),
"Oui"="1","Non"="2",NULL = "6666")
freq(coco.activite$risque.conflit.empl)Poursuivez les recodages…
# Recodages suivants :
coco.activite$risque.conflit.coll <- fct_recode(factor(coco.activite$q41_3),
"Oui"="1","Non"="2",NULL = "6666")
coco.activite$risque.cho <- fct_recode(factor(coco.activite$q41_4),
"Oui"="1","Non"="2",NULL = "6666")
coco.activite$risque.baisse.sal <- fct_recode(factor(coco.activite$q41_5),
"Oui"="1","Non"="2",NULL = "6666")
coco.activite$risque.aucun <- fct_recode(factor(coco.activite$q41_6),
"Oui"="1","Non"="2",NULL = "6666")Représentez une de ces variables en graphique :
coco.activite %>%
filter(is.na(risque.sante)==F) %>%
ggplot()+
aes(fill=risque.sante, x=sexe)+
geom_bar(position="fill")Ensuite ?
On essayera de présenter tous ces croisements “collés” dans un unique tableau.
rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.sante))## risque
## sexe Oui Non Total
## Homme 47.8 52.2 100.0
## Femme 63.7 36.3 100.0
## Ensemble 57.3 42.7 100.0
On utilise rbind qui “colle” 2 tables en lignes. Avec 2
risques :
risques <- rbind(
# 1er tableau croisé :
rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.sante)) %>%
data.frame() %>%
mutate(type="Santé"),
# 2ème tableau croisé :
rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.conflit.empl)) %>%
data.frame() %>%
mutate(type="Conflit employeur")) Puis mise en forme et graphique :
risques %>%
filter(risque=="Oui" & sexe !="Ensemble") %>%
ggplot()+
aes(x=type, y=Freq)+
geom_bar(stat="identity")+
facet_wrap(~sexe)Poursuivre avec l’ensemble des risques
risques <- rbind(risques,
rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.conflit.coll)) %>%
data.frame() %>%
mutate(type="Conflit collègues"))
risques <- rbind(risques,
rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.cho)) %>%
data.frame() %>%
mutate(type="Chômage"))
risques <- rbind(risques,
rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.baisse.sal)) %>%
data.frame() %>%
mutate(type="Baisse de salaire"))
risques <- rbind(risques,
rprop(table(sexe=coco.activite$sexe, risque=coco.activite$risque.aucun)) %>%
data.frame() %>%
mutate(type="Aucun"))Et le graphique :
risques %>%
filter(risque=="Oui" & sexe !="Ensemble") %>%
ggplot()+
aes(x=type, fill=sexe, y=Freq)+
geom_bar(stat="identity",position="dodge")+
theme(axis.text.x = element_text(angle=45, hjust=1))+
labs(title="Exposition aux risques professionnels", fill="Genre",
x="",y="%",
caption = "Champ : Personnes en activité\nSource : Coco1 - ELIPSS, CDSP, avril 2020")Et la table :
risques %>%
filter(risque =="Oui") %>%
mutate(Freq = round(Freq,1)) %>%
pivot_wider(-risque,names_from = "sexe", values_from ="Freq")On l’a également vu l’an dernier mais un rapide rappel de la façon de joindre différentes tables vous aidera pour le dossier d’évaluation.
L’enquête “Faire face au covid19” s’est faite en plusieurs vagues. Tou·tes les enquêté·es n’ont pas répondu à chaque vague, il faudra donner certaines précisions sur le type de jointure attendu.
coco2 <- read.csv("data/cdsp_coco2_2020/donnees_coco2/csv_coco2/fr.cdsp.ddi.elipss.202004.csv",
stringsAsFactors = F)
jointure.left <- left_join(coco, coco2,
by = "UID_ANONYM_COVID")
jointure.right <- right_join(coco, coco2,
by = "UID_ANONYM_COVID")
jointure.full <- full_join(coco, coco2,
by = "UID_ANONYM_COVID")
jointure.inner <- inner_join(coco, coco2,
by = "UID_ANONYM_COVID")
nrow(jointure.left)## [1] 1076
nrow(jointure.right)## [1] 998
nrow(jointure.full)## [1] 1166
nrow(jointure.inner)## [1] 908
D’ici la prochaine séance, si vous souhaitez poursuivre les exercices, vous pouvez reprendre le déroulé du corrigé de l’examen de l’an dernier (que je vous avais déjà envoyé l’an dernier en M1) -> Section “M1 R Initiation (2020/2021)” > Examen. On pourra y revenir rapidement en début de prochaine séance si vous avez des questions.
D’autres recodages avec {tidyverse} : https://juba.github.io/tidyverse/09-recodages.html
Les graphiques avec {ggplot2} : https://juba.github.io/tidyverse/08-ggplot2.html
# Chargement du package :
library(tidyverse)
# Importation des données :
etab <- read.csv2("data/Examen_M1_21-22/fr-esr-principaux-etablissements-enseignement-superieur.csv",
stringsAsFactors = F, fileEncoding = "UTF-8")
etu <- read.csv2(
"data/Examen_M1_21-22/fr-esr-statistiques-sur-les-effectifs-d-etudiants-inscrits-par-etablissement.csv",
stringsAsFactors = F, fileEncoding = "UTF-8")
perm <- read.csv2("data/Examen_M1_21-22/fr-esr-enseignants-titulaires-esr-public.csv",
stringsAsFactors = F, fileEncoding = "UTF-8")
nperm <- read.csv2("data/Examen_M1_21-22/fr-esr-enseignants-nonpermanents-esr-public.csv",
stringsAsFactors = F, fileEncoding = "UTF-8")La principale difficulté a été la non-prise en compte de l’encodage. !! Vérifiez toujours que les données ont correctement été importées. Les aspects à prendre en compte :
-> Conditionne la fonction à utiliser pour lire le contenu du fichier et le stocker dans un objet
-> Préciser si besoin le caractère séparateur et Vérifier avec
names(table) ou str(table)
-> En version numérique, les caractères sont stockés sous forme de codes en octets. Notamment en raison de différents alphabets, il existe plusieurs systèmes de reconnaissance de ces codes : C’est l’encodage. En plus de ça, les systèmes d’exploitation ont des systèmes d’encodage par défaut différents : Linux / Mac <-> UTF-8 ; Windows (ordinateur configuré en langue Française) <-> WINDOWS-1252 (ou ISO-8859-1). UTF-8 est plus universel et est plus standard. Beaucoup de données sont fournies avec ce système d’encodage.
-> Préciser le système d’encodage et vérifier que les caractères
spéciaux (accentués, par exemple) sont correctement lisibles avec
unique(table$VariableNominale) par exemple.
Table Établissements de l’enseignement supérieur en France :
Quelle était la difficulté ?
class(etab$Effectifs.d.étudiants.inscrits.2017.18)## [1] "character"
head(etab$Effectifs.d.étudiants.inscrits.2017.18)## [1] "2 333" "1 101" "" "986" "6 103" "2 154"
La colonne contenant le nombre d’étudiant·es dans les établissements n’était pas de type numérique. De plus, elle comportait des espaces insécables séparant les milliers. Il arrive que les données obtenues ne soient pas formatées comme on l’aurait souhaité. Comment la reformater correctement ?
etab$etu_1718 <- gsub(" ","",etab$Effectifs.d.étudiants.inscrits.2017.18)
etab$etu_1718 <- as.integer(etab$etu_1718)
mean(etab$etu_1718, na.rm=T)## [1] 8593.155
En version R-base :
tapply(etab$etu_1718, etab$secteur_d_etablissement, mean, na.rm=T)## Privé Public
## NaN 2686.532 12047.972
En version tidyverse :
etab %>%
group_by(Secteur=secteur_d_etablissement) %>%
summarise(`Nb moyen d'étudiant·es en 17/18` = mean(etu_1718, na.rm=T))On aurait aussi pu ajouter le nombre d’établissements pour chacun des secteurs et arrondir la moyenne :
etab %>%
group_by(Secteur=secteur_d_etablissement) %>%
summarise(`Nb moyen d'étudiant·es en 17/18` = round(mean(etu_1718, na.rm=T),0),
"Nombre d'établissements"=n())etab$date_creation_cl <- case_when(is.na(etab$date_creation) ~ "NR",
etab$date_creation < 1900 ~ "Avant 1900",
etab$date_creation < 1945 ~ "1900-1944",
etab$date_creation < 1970 ~ "1945-1969",
etab$date_creation < 2000 ~ "1970-1999",
T ~ "2000 ou après")
# Mettre « Avant 1900 » en 1er :
etab$date_creation_cl <- fct_relevel(etab$date_creation_cl, "Avant 1900")
library(questionr)
freq(etab$date_creation_cl)freq(etab$compte_youtube!="")addmargins(table(etab$Ancienne.région, etab$secteur_d_etablissement))##
## Privé Public Sum
## 0 0 6 6
## Alsace 0 0 7 7
## Aquitaine 0 2 6 8
## Auvergne 0 0 2 2
## Basse-Normandie 0 1 3 4
## Bourgogne 0 1 2 3
## Bretagne 0 5 13 18
## Centre-Val de Loire 0 0 3 3
## Champagne-Ardenne 0 3 3 6
## Corse 0 0 1 1
## Franche-Comté 0 1 3 4
## Guadeloupe 0 0 1 1
## Guyane 0 0 1 1
## Haute-Normandie 0 3 6 9
## Île-de-France 1 35 53 89
## La Réunion 0 0 1 1
## Languedoc-Roussillon 0 1 4 5
## Limousin 0 1 1 2
## Lorraine 0 2 1 3
## Mayotte 0 0 1 1
## Midi-Pyrénées 0 3 16 19
## Nord-Pas-de-Calais 0 7 5 12
## Pays de la Loire 0 11 4 15
## Picardie 0 4 3 7
## Poitou-Charentes 0 2 3 5
## Provence-Alpes-Côte d'Azur 0 2 7 9
## Rhône-Alpes 0 5 14 19
## Sum 1 89 170 260
rprop(table(etab$Ancienne.région, etab$secteur_d_etablissement))##
## Privé Public Total
## 0.0 0.0 100.0 100.0
## Alsace 0.0 0.0 100.0 100.0
## Aquitaine 0.0 25.0 75.0 100.0
## Auvergne 0.0 0.0 100.0 100.0
## Basse-Normandie 0.0 25.0 75.0 100.0
## Bourgogne 0.0 33.3 66.7 100.0
## Bretagne 0.0 27.8 72.2 100.0
## Centre-Val de Loire 0.0 0.0 100.0 100.0
## Champagne-Ardenne 0.0 50.0 50.0 100.0
## Corse 0.0 0.0 100.0 100.0
## Franche-Comté 0.0 25.0 75.0 100.0
## Guadeloupe 0.0 0.0 100.0 100.0
## Guyane 0.0 0.0 100.0 100.0
## Haute-Normandie 0.0 33.3 66.7 100.0
## Île-de-France 1.1 39.3 59.6 100.0
## La Réunion 0.0 0.0 100.0 100.0
## Languedoc-Roussillon 0.0 20.0 80.0 100.0
## Limousin 0.0 50.0 50.0 100.0
## Lorraine 0.0 66.7 33.3 100.0
## Mayotte 0.0 0.0 100.0 100.0
## Midi-Pyrénées 0.0 15.8 84.2 100.0
## Nord-Pas-de-Calais 0.0 58.3 41.7 100.0
## Pays de la Loire 0.0 73.3 26.7 100.0
## Picardie 0.0 57.1 42.9 100.0
## Poitou-Charentes 0.0 40.0 60.0 100.0
## Provence-Alpes-Côte d'Azur 0.0 22.2 77.8 100.0
## Rhône-Alpes 0.0 26.3 73.7 100.0
## Ensemble 0.4 34.2 65.4 100.0
# En nombres en Île-de-France – où se trouvent 1/3 des établissements
# de France (freq(etab$Ancienne.région=="Île-de-France")) -, en proportion
# dans les Pays de la Loire.etab %>%
filter(Ancienne.région =="Nord-Pas-de-Calais") %>%
arrange(desc(etu_1718)) %>%
select(Nom=Libellé, Secteur=secteur_d_etablissement,
`Nb d'étudiant·es`=etu_1718, `Date de création`=date_creation) %>%
slice(1:5)Table des Effectifs étudiants
etu <- etu %>% filter(Année.universitaire=="2020-21")# names(etu)
etu <- etu %>%
rename(etu_total =
Nombre.d.étudiants.inscrits..inscriptions.principales..y.compris.doubles.inscriptions.CPGE) %>%
mutate(ratioFemmes = round(Sexe...Femmes/etu_total*100,2),
ratioInge =
round(Grande.discipline...Sciences.et.sciences.de.l.ingénieur/etu_total*100,2))etu %>%
filter(ratioInge >0) %>%
ggplot() +
aes(x = ratioFemmes, y = ratioInge, size = etu_total,
col = Type.d.établissement) +
geom_point(alpha=.6) +
theme_bw()+
labs(title = "Femmes en sciences de l'ingénierie - 20/21",
x = "Part de femmes (%)", y = "Part en sciences de l'ingénierie (%)",
caption = "Source : Ministère de l'ESRI, 2021",
size = "Nombre d'étudiant·es",
col = "Type d'étab.")full_join(perm %>%
filter(Rentrée == "2019") %>%
filter(Région=="Hauts-de-France") %>%
group_by( Établissement) %>%
summarise(permanents=sum(effectif)),
nperm %>%
filter(Rentrée == "2019") %>%
filter(Région=="Hauts-de-France") %>%
group_by( Établissement) %>%
summarise(`non-permanents`=sum(Effectif)),
by="Établissement") %>%
mutate(`Part de non-permanent·es` =
round(`non-permanents`/(`non-permanents`+permanents)*100,2)) %>%
arrange(desc(`Part de non-permanent·es`)) Source : allison
Horst @allison_horst
Les données sont organisées sous forme de couches superposables. On distingue généralement deux types de données : les données de type vecteur et les raster.
Les données vecteurs se définissent uniquement par des coordonnées. On trouvera des données vecteurs de type points, lignes et polygones. Un point sera défini par un couple de coordonnées XY, une ligne ou un polygone par les coordonnées de leurs sommets. Une couche vecteur sera soit de type point, soit de type ligne, soit de type polygone. Des données attributaires - supplémentaires - sont associées aux données spatiales : le nombre d’habitants d’une ville, par exemple.
Les données raster sont des images (vue satellite, carte administrative, carte IGN).
La Terre est (à peu près) ronde et se présente en 3 dimensions. Afin de représenter les informations géographiques sur une carte en 2 dimensions, il faudra projeter les coordonnées spatiales, ce qui donnera lieu à des déformations. Les multiples méthodes de projections ont leurs avantages et inconvénients et seront utilisées selon différentes finalités.
Par exemple, la projection de Mercator conserver les angles mais pas les distances. Les surfaces et distances s’éloignant de l’équateur paraissent plus importantes qu’elles ne le sont en réalité. En revanche, cette projection conserve les angles et se trouve donc utile en navigation.
Extraits de “Le scandale de Mercator”, les Savoirs ambulants, https://www.lessavoirsambulants.fr/p/le-scandale-de-mercator.html
Un type de projection nous intéressera, il s’agit des projections coniques qui minimisent les déformations sur des zones du globes.
Utilisée en cartographie, on représente différentes zones en modifiant le référentiel. Le référentiel Lambert-93 est adapté pour représenter la France métropolitaine. Les systèmes de coordonnées sont multiples. Vous en rencontrerez deux principaux : WGS84 et Lambert-93 identifiables par les EPSG (European Petroleum Survey Group) 4326 et 2154.
Il est essentiel de s’assurer que tous vos fichiers spatiaux appartiennent au même système de coordonnées de référence (SCR). Les fichiers spatiaux peuvent être re-projetés automatiquement.
sf est un des packages de référence pour traiter
simplement des données spatiales dans R. Il a l’avantage d’être
compatible avec l’esprit et les opérateurs tidyverse.
/!\ Comme tout le reste dans R, les choses évoluent, et rapidement. Notamment pour les packages les plus utilisés et maintenus ainsi que les formats de données les plus répandus.
Dès maintenant, installez les packages suivants (qui sont assez longs à installer, ne laissez pas cette ligne directement dans le script) :
install.packages(c("sf", "geojsonio","ggspatial"))Puis, téléchargez le fichier “carto.zip” du moodle et mettez-le dans le dossier data :
Exemple avec des premiers jeux de données sur la page d’exploration des données open data de la MEL :
Créer un dossier et extraire le contenu d’une archive zip :
# Créer un dossier "carto" dans le dossier data s'il n'existe pas déjà :
if(dir.exists("data/carto")==F){
dir.create("data/carto")
}
# Extraire les fichiers du zip dans le dossier créé :
unzip(zipfile="data/carto.zip", exdir = "data/carto")Shapefile = Le format shapefille est une des formes de stockage de données vectorielles qui s’est imposée (utilisée dans le SIG ArcGIS et compatible avec des logiciels libres comme QGIS). Il est en fait composé de plusieurs fichiers :
library(sf)
# Import du fichier shp :
garesMEL_sf <- st_read(dsn="data/carto/liste-des-gares.shp")## Reading layer `liste-des-gares' from data source
## `/home/cecile/Documents/Cours et ateliers/ENSP_22-23/M2 - Approfondissement R/data/carto/liste-des-gares.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 59 features and 15 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 2.804054 ymin: 50.52747 xmax: 3.234699 ymax: 50.76164
## Geodetic CRS: WGS 84
class(garesMEL_sf)## [1] "sf" "data.frame"
Vous trouverez aussi parfois des fichiers vectoriels en format GeoJSON. Pour les lire :
library(geojsonio)
# Importer un fichier GeoJSON :
communesMEL_sf <- geojson_read("data/carto/mel_communes.geojson", what="sp") %>%
st_as_sf() # Le transformer en format sfLes objets spatiaux sf sont constitués d’une partie
“géométrique” (“geometry”, informations spatiales) et d’une partie
“données” (attributs qui seront représentés) sous la forme d’un
data.frame.
names(communesMEL_sf)## [1] "ut" "nom" "insee" "territoire" "surface"
## [6] "code_posta" "perimetre" "geo_point_2d" "geometry"
names(garesMEL_sf)## [1] "code_uic" "libelle" "fret" "voyageurs" "code_ligne"
## [6] "rg_troncon" "pk" "commune" "departemen" "idreseau"
## [11] "idgaia" "x_l93" "y_l93" "x_wgs84" "y_wgs84"
## [16] "geometry"
La partie geometry peut être représentée comme un type de graphique :
# L'objet `data` est précisé dans chaque geom_sf ou couche cartographique
ggplot()+
geom_sf(data=communesMEL_sf)Il est important de vérifier le sytèmes de coordonnées : les graphiques reprojettent les couches à la volée mais certaines opérations ne seront pas possibles s’ils ne sont pas compatibles.
# st_crs(communesMEL_sf)
# st_crs(garesMEL_sf)
communesMEL_sf <- st_transform(communesMEL_sf,crs=2154)
garesMEL_sf <- st_transform(garesMEL_sf,crs=2154)
# st_crs(communesMEL_sf)
# st_crs(garesMEL_sf)Représenter plusieurs couches sur une carte :
ggplot()+
geom_sf(data=communesMEL_sf)+
geom_sf(data=garesMEL_sf)La carte se construit donc par ajout d’éléments esthétiques, comme pour une graphique de données non-géographiques :
ggplot()+
geom_sf(data=communesMEL_sf, aes(fill=surface))+
viridis::scale_fill_viridis()+
geom_sf(data=garesMEL_sf,aes(col=voyageurs))+
scale_color_manual(values = c("black","grey"))+
theme_bw()+
labs(title = "Surface des communes et gares de la MEL", col="Voyageurs",
fill="Surface (m²)", caption="Source : MEL")+
coord_sf( datum = NA)Les données directement incluses dans les fichiers spatiaux ne sont pas toujours suffisantes et gagnet à être enrichies. Les éléments (points, polygones, lignes) possèdent souvent des identifiants. Par exemple, pour le fichier des communes de la MEL, la colonne “insee” correspond aux codes INSEE unique des communes (attentions, il ne s’agit pas des codes postaux).
Importer les données de la population légale
des communes du recensement 2018 (fichier Communes.csv) et ajouter les
informations au fichier spatial à l’aide d’une jointure de table avec
left_join.
pop18 <- read.csv2("data/carto/Communes.csv") # import
communesMEL_sf <- left_join(communesMEL_sf, pop18, by =c("insee" = "DEPCOM")) # jointure
communesMEL_sf <- communesMEL_sf %>% filter(is.na(PTOT)==FALSE)Représenter la population de chaque commune.
ggplot()+
geom_sf(data=communesMEL_sf,aes(fill=PTOT),color="snow4",size=.2)+
viridis::scale_fill_viridis()+
theme_bw()+
labs(title = "Population des communes de la MEL",
fill="Population", caption="Source : RP 2017, Insee")+
coord_sf( datum = NA)Représenter les communes selon les quintiles
de population. On crée une variable représentant les quintiles de
population avec les fonctions cut et
quantile
communesMEL_sf$PTOT_quintiles <- cut(communesMEL_sf$PTOT,
quantile(communesMEL_sf$PTOT, seq(0,1,.2), na.rm=T))
ggplot()+
geom_sf(data=communesMEL_sf,aes(fill=PTOT_quintiles),
color="snow4",size=.2)+
viridis::scale_fill_viridis(discrete=T)+
theme_bw()+
labs(title = "Population des communes de la MEL",
fill="Population", caption="Source : RP 2017, Insee")+
coord_sf( datum = NA)Au-delà de données contenues dans les fichiers spatiaux et de données extérieures les enrichissant, toute une série de transformation et de calculs peuvent être directement effectués sur les données spatiales : Combien y’a-t-il de gares dans chaque commune de la MEL ? Quel est le prix moyen d’achat des logements dans chaque quartier de Lille ?
Les fonctions de transformation et de calculs sur les fichiers spatiaux apparaissent avec le préfixe “st_”.
Les données DVF sont ” publiées et produites par la direction générale des finances publiques, permettent de connaître les transactions immobilières intervenues au cours des cinq dernières années sur le territoire métropolitain et les DOM-TOM, à l’exception de l’Alsace-Moselle et de Mayotte. Les données contenues sont issues des actes notariés et des informations cadastrales.” (https://www.data.gouv.fr/fr/datasets/5c4ae55a634f4117716d5656/)
Importer les données DVF (59.csv)
dvf_brut <- read.csv("data/carto/59.csv", stringsAsFactors = F)Les données DVF sont des données issues des actes notariaux et pas forcément utilisables directement. Il faut passer par des traitements de nettoyage. Pour gagner du temps, j’ai préparé les données et fait les modifications suivantes :
Le fichier nettoyé à importer est “dvf19_cours.rds”.
dvf <- readRDS("data/carto/dvf19_cours.rds")Pour autant, ce n’est pas un fichier spatial … Mais on dispose de 2 colonnes contenant les informations “latitude” et “longitude”. On peut s’en servir pour créer un fichier spatial. On construit la partie “géométrique” à l’aide des colonnes longitude et latitude et de la fonction st_as_sf(). Les coordonnées sont exprimées selon la projection WGS84 (crs = 4326)
On peut transformer ce tableau de données en un fichier spatial à
l’aide de st_as_sf en indiquant les colonnes de longitude
(x) et latitude (y) et le crs (ici, 4326) :
dvf_sf <- dvf %>%
st_as_sf(coords = c("longitude", "latitude"), crs=4326)
dvf_sf <- dvf_sf %>% st_transform(crs=2154)Puis, on représente les points des ventes et les délimitations des communes de la MEL :
ggplot() +
geom_sf(data=dvf_sf) +
geom_sf(data=communesMEL_sf)Pour la suite, on s’intéresse aux ventes immobilières à une échelle plus fine, celle des IRIS (unité territoriale infracommunale pour les communes de plus de 10.000 habitants) de la MEL (Métropole Européenne de Lille, fichier MEL_iris.shp). Importer le fichier spatial “MEL_iris” puis représentez les ventes immobilières à l’échelle des IRIS sur une carte.
irisMEL_sf <- st_read("data/carto/MEL_iris.shp")## Reading layer `MEL_iris' from data source
## `/home/cecile/Documents/Cours et ateliers/ENSP_22-23/M2 - Approfondissement R/data/carto/MEL_iris.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1271 features and 18 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 685041 ymin: 7044709 xmax: 719291.4 ymax: 7077561
## Projected CRS: RGF93_Lambert_93
# st_crs(MEL_iris)
ggplot() +
geom_sf(data=dvf_sf) +
geom_sf(data=irisMEL_sf)Pour associer chaque point à chaque iris, la fonction
st_intersection permet de sortir la zone commune à deux
fichiers sf en ajouter les données des iris au fichier des points des
ventes immobilières.
dvf_MEL_sf <- st_intersection( dvf_sf, irisMEL_sf)ggplot() +
geom_sf(data=dvf_MEL_sf)+
geom_sf(data=irisMEL_sf, fill="transparent")On a oublié de regarder les systèmes de projection de ces 2 objets ! Ayez ce réflexe, même si les fonctions de carto reprojetent les données à la volée.
Un fichier sf se compose d’un tableau de données et des informations spatiales associées à chaque entité.
Produire une table contenant le nombre de
transactions immobiliers par IRIS (CODE_IR) à l’aide des fonctions
group_by et summarise.
transactions_IRIS <- dvf_MEL_sf %>%
group_by(CODE_IR) %>%
summarise(nb_transactions=n()) %>%
data.frame() %>%
select(CODE_IR, nb_transactions)Représenter cette données sur une carte (le fonds de carte des IRIS est irisMEL_sf)
irisMEL_sf <- left_join(irisMEL_sf, transactions_IRIS, by="CODE_IR")
ggplot()+ geom_sf(data=irisMEL_sf,size=.2,
aes(fill=nb_transactions))+
viridis::scale_fill_viridis()+
theme_bw()+
coord_sf(crs = st_crs(irisMEL_sf), datum = NA)Représenter la taille moyenne (variable surface_reelle_bati) des logements vendus par IRIS dans la MEL puis à Lille.
taille <- dvf_MEL_sf %>%
# filter(surface_reelle_bati<800) %>%
group_by(CODE_IR) %>%
summarise(taille = mean(surface_reelle_bati)) %>%
data.frame() %>%
select(CODE_IR, taille)
irisMEL_sf$taille <- NULL
irisMEL_sf <- left_join(irisMEL_sf, taille, by="CODE_IR")
ggplot()+
geom_sf(data=irisMEL_sf, size=.1)+
aes(fill=taille)+
viridis::scale_fill_viridis()+
theme_bw()+
labs(title = "Taille des logements vendus dans la métropole lilloise",
fill="Surface réelle du bâti (m²)",
caption="Source : DVF, Ministère des Finances, 2019")irisLille_sf <- irisMEL_sf %>% filter(NOM_COM=="Lille")
ggplot()+
geom_sf(data=irisLille_sf, size=.1)+
aes(fill=taille)+
viridis::scale_fill_viridis()+
theme_bw()+
labs(title = "Taille des logements vendus à Lille",
fill="Surface réelle du bâti (m²)",
caption="Source : DVF, Ministère des Finances, 2019")Représenter les lieux des ventes dans Lille et les prix au m² de ces ventes.
dvf_Lille_sf <- st_intersection( dvf_MEL_sf, irisLille_sf)
dvf_Lille_sf$prix <- round(dvf_Lille_sf$valeur_fonciere/dvf_Lille_sf$surface_reelle_bati)
ggplot()+
geom_sf(data=irisLille_sf, size=.1, fill="transparent")+
geom_sf(data=dvf_Lille_sf, aes(col=prix), size=.1)+
viridis::scale_color_viridis()+
coord_sf(crs = st_crs(irisLille_sf), datum = NA)+
theme_bw()Le package ggspatial permet d’ajouter un fond de type
raster (image ou “tuiles”).
library(ggspatial)
ggplot()+
annotation_map_tile(zoom=13,type="cartolight")+
geom_sf(data=irisLille_sf, size=.3, fill="transparent", col="grey60")+
geom_sf(data=dvf_Lille_sf, aes(col=prix), size=.1)+
viridis::scale_color_viridis()+
theme_bw()+
labs(title = "Prix au m² des ventes de 2019 à Lille",
col="Prix (/m²)",
caption="Sources : DVF, Ministère des Finances, 2019, fond : OSM")+
coord_sf(crs = st_crs(irisLille_sf), datum = NA)Si on souhaite faire apparaître certains éléments de fond au lieu des tuiles images, on peut faire une extraction depuis OpenStreetMap à l’aide du package {osmdata} :
library(osmdata)
roads <- opq(bbox = st_bbox(irisLille_sf) ) %>%
add_osm_feature(key = 'highway') %>%
osmdata_sf()
roads <- roads$osm_lines
st_crs(roads) <- 4326
roadsLille_sf <- st_intersection( roads, irisLille_sf)
ggplot()+geom_sf(data=roads)
ggplot()+geom_sf(data=roadsLille_sf)Plus de détail sur le package {osmdata} : https://docs.ropensci.org/osmdata/
Ensemble des clés et valeurs d’interrogation de l’API Overpass turbo : https://wiki.openstreetmap.org/wiki/Map_features
Tutoriels sur les données spatiales avec R : https://rspatialdata.github.io/osm.html
Enfin, les cartes peuvent être rendues interactives ; C’est notamment utile si vous en fait une application {shiny} (voir séminaire/TD dédié). Plusieurs packages le permettent, vous pouvez vous reporter à ces deux articles :
Un exemple avec {tmap} :
library(tmap)
library(viridis)
tmap_mode(mode = "view") # Active le mode interactif
tm_basemap("CartoDB.Positron")+ # tuiles rasters
tm_shape(dvf_Lille_sf)+ # données
tm_symbols(col = "prix", border.col = "transparent",
scale = .3, alpha = 1, palette =viridis(n=10) ,
perceptual = TRUE) # Taille des points reste la mêmeReprésenter de la façon la plus lisible le prix au m² selon les catégories de taille des ventes de la MEL. Essayez de trouver une représentation graphique originale et adaptée à ce que vous voulez raconter !
# Prix au m² :
dvf_MEL_sf$prix <- round(dvf_MEL_sf$valeur_fonciere/dvf_MEL_sf$surface_reelle_bati)
# Catégories de tailles :
dvf_MEL_sf$taille <- case_when(dvf_MEL_sf$surface_reelle_bati<40 ~ "40-",
dvf_MEL_sf$surface_reelle_bati<=60 ~ "40-60",
dvf_MEL_sf$surface_reelle_bati<=80 ~ "60-80",
dvf_MEL_sf$surface_reelle_bati<=100 ~ "80-100",
dvf_MEL_sf$surface_reelle_bati<=120 ~ "100-120",
dvf_MEL_sf$surface_reelle_bati>120 ~ "120+",
T~"")
dvf_MEL_sf$taille <- factor(dvf_MEL_sf$taille,
levels=c("40-","40-60","60-80","80-100",
"100-120","120+"))# Graphique 1 :
# TaillePrix
TaillePrix <- dvf_MEL_sf %>%
group_by(taille) %>%
summarise("Prix moyen"=mean(prix),
"Prix médian"=median(prix))
TaillePrix %>%
ggplot()+
geom_linerange(aes(ymax = `Prix moyen`,
ymin=`Prix médian`,
x=taille), col="grey60", size=1.5)+
geom_point(aes(y= `Prix moyen`, x=taille), col="white",
size=1.5)+
geom_point(aes(y= `Prix moyen`, x=taille), col="firebrick",
shape=1, stroke=2, size=1.5)+
geom_point(aes(y= `Prix médian`, x=taille), col="white",
size=1.5)+
geom_point(aes(y= `Prix médian`, x=taille), col="seagreen",
shape=1, stroke=2, size=1.5)+
theme_bw()+ coord_flip()+
labs(y="Prix (rouge=moyen, vert=médian)",
x="Taille (m²)",
title="Prix/m² des logements selon la taille",
subtitle="Métropole Européenne de Lille (2019)",
caption="Source : DVF")# Graphique 2 :
library(ggridges)
dvf_MEL_sf %>% ggplot(aes(x=prix, y=taille, fill=stat(x)))+
geom_density_ridges_gradient(col="white",scale = 3, size = .5, rel_min_height = 0.01) +
scale_fill_viridis_c(name = "", option = "C")+
theme_bw()# https://www.datanovia.com/en/blog/elegant-visualization-of-density-distribution-in-r-using-ridgeline/Retour aux données Coco sur lesquelles vous allez travailler pour le rendu de fin de semestre :
Il se trouve que les données Coco sont des données construites sur échantillon. On l’a déjà entre-aperçu, la structure de la population est particulière alors que le panel se veut représentatif de la population de France métropolitaine. Pour corriger les effets de structure de l’échantillon enquêté, les producteur·rices des données ont construit des coefficients de redressement que nous allons utiliser.
On a vu que le panel était plutôt âgé. On retrouve cette correction dans les poids :
coco %>%
group_by(age_quinq) %>%
summarise(`poids médian`=median(POIDS),
`poids moyen`=mean(POIDS))Les plus jeunes sont “sur-pondéré·es” et les plus âgé·es “sous-pondéré·es”. Et la somme des POIDS est égale à la taille de l’échantillon, on dit que la pondération est normée.
Plusieurs packages permettent d’utiliser les pondérations :
wtd.table du package
{questionr}.{survey} est plus délicat à prendre en main que d’autres
packages mais il intègre l’ensemble des traitements et tests qu’on peut
vouloir faire. Il faut parfois un peu “bidouiller” des tableaux de
résultats pour les présenter.On commence par faire le tri à plat redressé des PCS. Il faut déjà recoder les PCS de la variable coco1_q35 :
coco$pcs <- fct_collapse(factor(coco$q35),
"Agriculteurs"="1",
"Artisans, comm., chef d'E"="2",
"Cadres et pr.intell.sup."="3",
"Prof. intermédiaires"="5",
"Employés"="4","Ouvriers"="6",
NULL = c("7","9999"))
freq(coco$pcs)# N pondéré :
wtd.table(coco$pcs, weights = coco$POIDS,
useNA = "ifany")## Agriculteurs Artisans, comm., chef d'E Cadres et pr.intell.sup.
## 26.78782 82.68868 222.98026
## Employés Prof. intermédiaires Ouvriers
## 167.20659 448.01096 87.77051
## <NA>
## 40.55518
Du fait des poids, les effectifs ont des décimales donc on les arrondit :
wtd.table(coco$pcs, weights = coco$POIDS,
useNA = "ifany") %>%
round()#<<## Agriculteurs Artisans, comm., chef d'E Cadres et pr.intell.sup.
## 27 83 223
## Employés Prof. intermédiaires Ouvriers
## 167 448 88
## <NA>
## 41
# % pondérés :
round(wtd.table(coco$pcs, weights = coco$POIDS,
useNA = "ifany") %>%
prop.table()*100,1)#<<## Agriculteurs Artisans, comm., chef d'E Cadres et pr.intell.sup.
## 2.5 7.7 20.7
## Employés Prof. intermédiaires Ouvriers
## 15.5 41.6 8.2
## <NA>
## 3.8
Croisons la PCS avec la variable q38 (Situation actuelle vis-à-vis du télétravail) :
# Recodages :
coco$q38 <- fct_recode(factor(coco$q38),
"Lieu habituel"="1",
"Alternance lieu habituel / TT"="2",
"Télé-travail"="3", NULL="6666")
# Tableau :
wtd.table(coco$pcs, coco$q38,
weights = coco$POIDS, useNA = "ifany") %>%
rprop(digits=1)## Lieu habituel Alternance lieu habituel / TT
## Agriculteurs 42.2 3.0
## Artisans, comm., chef d'E 2.8 1.0
## Cadres et pr.intell.sup. 4.4 3.3
## Employés 14.0 2.6
## Prof. intermédiaires 18.5 1.5
## Ouvriers 12.0 0.0
## <NA> 0.0 0.0
## Ensemble 13.1 1.8
## Télé-travail <NA> Total
## Agriculteurs 0.0 54.7 100.0
## Artisans, comm., chef d'E 4.5 91.6 100.0
## Cadres et pr.intell.sup. 27.6 64.6 100.0
## Employés 12.2 71.2 100.0
## Prof. intermédiaires 5.6 74.4 100.0
## Ouvriers 0.8 87.2 100.0
## <NA> 10.9 89.1 100.0
## Ensemble 10.8 74.3 100.0
{lourdR} est un package que j’ai commencé à faire. Il est perfectible mais utile. Je l’ai déposé sur github, pour l’installer, voici les lignes à exécuter (une fois) :
if (!require("devtools")) install.packages("devtools", dep=T)
# installe devtools si vous ne l'avez pas
devtools::install_github("Grisoudre/lourdR") # installe lourdREnsuite, il peut être chargé comme un package classique :
library(lourdR){lourdR} contient - pour le moment - deux fonctions, on utilise déjà
freqp : copie de la fonction
freq présentant les résultats non-pondérés et pondérés
:
freqp(coco, var = "pcs", poids = "POIDS")# n=F : sans les effectifs
# p=F : sans les pourcentages
# brut=F : sans les valeurs brutes
# exclude = "item" : Valeur à exclureOn peut voir l’effet des redressements sur les résultats !
L’autre fonction freqm présente dans un
même tableau des colonnes avec un même préfixe et des mêmes modalités -
typique des variables issues de questions à choix
multiples. Prenons les variables q24_1 à
q24_7 (Vit avec : …) :
freqm(coco, "q24", calc="%")
# calc = "n", "%", ou "all"
# exclude = "item" : valeur exclue
# poids = "variable_poids" : si données pondéréesPuis on renomme les modalités de “name”. On peut également pondérer ce tableau.
En retirant les “non-concerné·es”, en pondérant :
freqm(coco, "q24", calc="%", exclude="6666", poids="POIDS")Un peu mieux mis en forme avec {tidyverse} :
freqm(coco, "q24", calc="%", exclude="6666", poids="POIDS") %>%
mutate(name= fct_recode(name,
"Conjoint·e"="q24_1",
"Enfants/petits-enfants"="q24_2",
"Parents/Grands-parents"="q24_3",
"Famille (autre)"="q24_4",
"Ami·es ou autres"="q24_5",
"Chien(s)"="q24_6",
"Autre(s) animal(aux)"="q24_7")) %>%
rename("Vit avec"="name", "Oui"="X1") %>%
select(-X2, -Total) %>%
arrange(desc(Oui))Pour le moment, on a utilisé les poids pour présenter des résultats de types effectifs et pourcentages sur effectifs. C’est en fait assez limité. Si on souhaite faire des calculs plus avancés - types régressions - on peut se référer à un autre package, {survey}. Il est en général présenté dans les tutoriels et enseignements de R en SHS (https://larmarange.github.io/analyse-R/donnees-ponderees.html). Il est plus délicat à prendre en main que d’autres packages mais il intègre l’ensemble des traitements et tests qu’on peut vouloir faire. Il faut parfois un peu “bidouiller” des tableaux de résultats pour les présenter. Installons-le :
install.packages("survey", dep = T)library(survey)Avant toute chose, il faut créer un nouvel objet de type
“survey design” explicitant le plan d’échantillonnage
(ou design) et donc la variable de pondération à considérer. Donc on
produit le survey design avec svydesign
:
coco$sexe <- fct_recode(factor(coco$eayy_A1),"Homme"="1","Femme"="2")
# Recodage du genre avant
cocow <- svydesign(data = coco,
weights = ~ coco$POIDS,
ids = ~1)
class(cocow) # = Objet bizarre## [1] "survey.design2" "survey.design"
L’objet cocow n’est plus un tableau de données. Le
tableau de données se trouve en fait dans la partie “variables” du
“survey.design” :
names(cocow$variables)[1:10]## [1] "UID_ANONYM_COVID" "q01" "q02" "q03"
## [5] "q04" "q05_c" "q06_1" "q06_2"
## [9] "q06_3" "q06_4"
Donc !!!!! Il est important de faire tous vos recodages et sélections d’individus AVANT de créer l’objet survey design. N’oubliez pas qu’un script est linéaire.
On peut toujours faire des calculs sur le tableau initial :
freq(cocow$variables$sexe)On peut aussi être tenté·e de faire les recodages et modifications sur ce tableau-là (je le fais aussi) mais il faut être sûr·e de soi, ça complique vite les choses. Conseil : Faire les recodages, créations de variables, etc. avant !
svytable(~variable, design)svytable(~sexe, cocow)## sexe
## Homme Femme
## 529.2982 546.7018
–
svytable(~sexe, cocow, round = T)## sexe
## Homme Femme
## 529 547
freq) :svytable(~sexe, cocow, round = T) %>%
freq(valid=F, total = TRUE)svytable(~pcs, cocow, round=T) %>%
freq(valid=F, total = TRUE){survey} ignore les NA !!
–
Il faut les expliciter
cocow$variables$pcs <- fct_explicit_na(cocow$variables$pcs,"NR")
svytable(~pcs, cocow, round=T) %>%
freq(valid=F, total = TRUE)svytable(~pcs+q38, cocow) %>%
rprop(digits=1)## q38
## pcs Lieu habituel Alternance lieu habituel / TT
## Agriculteurs 93.3 6.7
## Artisans, comm., chef d'E 33.7 12.2
## Cadres et pr.intell.sup. 12.6 9.3
## Employés 48.6 8.9
## Prof. intermédiaires 72.3 5.8
## Ouvriers 93.5 0.0
## NR 0.0 0.0
## Ensemble 50.8 7.2
## q38
## pcs Télé-travail Total
## Agriculteurs 0.0 100.0
## Artisans, comm., chef d'E 54.1 100.0
## Cadres et pr.intell.sup. 78.2 100.0
## Employés 42.4 100.0
## Prof. intermédiaires 21.9 100.0
## Ouvriers 6.5 100.0
## NR 100.0 100.0
## Ensemble 42.0 100.0
Exemple d’une moyenne pondérée (q12 = nombre de personnes dans l’entourage ayant contracté le covid) :
svymean(~ q12, cocow, na.rm=T)## mean SE
## q12 2.8041 0.1465
{ggplot2} permet d’intégrer les données de survey design. Il faut ajouter deux précision :
ggplot()aes()C’est un peu du bidouillage…
Ex : Graphique du lieu d’activité professionnelle pendant le 1er confinement selon la PCS :
# Graphique :
ggplot(cocow$variables) +#<<
aes(weight = weights(cocow),#<<
x = pcs, fill = q38) +
geom_bar(position = "fill")+
theme(axis.text.x = element_text(angle=45, hjust=1))+
labs(title = "Lieu de travail au cours des deux dernières semaines",
x="PCS",fill="",y="",
caption ="Source : Coco1 - ELIPSS, CDSP, avril 2020")Exercices supplémentaires facultatifs : Refaire toutes les questions de la partie I. en version pondérée
Aller plus loin pour traiter les données pondérées avec {survey} : https://larmarange.github.io/analyse-R/donnees-ponderees.html
Nous verrons la façon de produire des tests statistiques sur des données pondérées dans la partie suivante.
Echantillon :
= Une partie de la population qu’on souhaite interroger. Les personnes interrogées – l’échantillon – sont choisies de façon à être représentatives de l’ensemble de la population : La manière dont sont choisies ces personnes est appelée échantillonnage.
Représentativité :
La représentativité statistique permettra de généraliser des résultats et observations faites sur l’échantillon à l’ensemble de la population. Pour être représentatif, un échantillon doit présenter les mêmes caractéristiques que sa population (soit par règle statistique d’aléatoire, soit par construction “raisonnée” d’un échantillon).
La représentativité statistique s’appuie sur les notions de probabilité. Elle permettra d’adopter un raisonnement déductif : C’est à dire de partir du particulier (l’échantillon) pour généraliser des résultats observés à un ensemble plus vaste (la population) = Statistiques inférentielles.
La statistique inférentielle consiste à généraliser des résultats observés sur un échantillon à l’ensemble de la population de référence. Produire des statistiques inférentielles se fait donc sous certaines conditions. Une des premières conditions nécessaires à la réalisation de statistiques inférentielles porte sur l’échantillon ; celui-ci doit pouvoir garantir une certaine représentativité statistique :
Avant de traiter les tests statistiques, voici deux exemples d’estimation de valeurs. Dans cette partie, à partir de données fictives, on découvrira comment estimer une valeur dans la population générale en utilisant qu’une partie des individus de cette population (un échantillon).
Exemples pour observer comment le choix de l’échantillon agit sur les résultats :
Soit, pour l’exemple 1, une population connue fictive de 10 000 individus qui ont mangé du chocolat l’année dernière.
On construit (fictivement) un tableau reprenant la quantité (en kg) de chocolat mangée l’an dernier par chacune de ces 10 000 personnes au cours de l’année écoulée. toutes ces données sont habituellement inconnues (c’est la quantité de chocolat consommée par toute la population d’un pays, par exemple, qu’on ne peut pas interroger de façon exhaustive).
dt <- data.frame(choco=round(rnorm(10000, mean = 13.2 , sd = 5)))
dt[dt<0]<-0
# loi normale de moyenne 13,2 et d'écart type 5
head(dt)L’individu avec l’identifiant 1 a mangé 21 kg de chocolat, la personne avec l’identifiant 2 en a mangé 11 kg.
La moyenne (habituellement inconnue et qu’on veut estimer) de cette quantité de chocolat ingurgitée est connue et de :
mean(dt$choco)## [1] 13.1801
= Valeur qu’on veut estimer !
Et les individus se répartissent comme ceci selon la quantité de chocolat mangée pendant l’année :
ggplot(dt)+
aes(x=choco) +
geom_histogram(stat="count")On reconnaît la forme de la répartition normale.
On demande à 500 d’entre elles et eux, l’échantillon, choisi·es au hasard (aléatoire simple), la quantité de chocolat mangée au cours des 12 derniers mois. On établit ensuite la moyenne de chocolat consommé par ces 500 individus, soit la moyenne estimée.
dt$sample <- sample(1:10000,10000,replace=F)
tirage <- data.frame(tirage=1,
moyenne=mean(dt[dt$sample <501,"choco"]),
sd = sd(dt[dt$sample <501,"choco"]),
n=nrow(dt[dt$sample <501,]))
tirageSur ce 1er tirage, on a une moyenne de 13.27.
On répète une 2ème fois l’opération avec 500 individus choisis au hasard parmi la population de départ ; on obtient une moyenne de :
dt$sample <- sample(1:10000,10000,replace=F)
tirage2 <- data.frame(tirage=2,
moyenne=mean(dt[dt$sample <501,"choco"]),
sd = sd(dt[dt$sample <501,"choco"]),
n=nrow(dt[dt$sample <501,]))
tirage <- rbind(tirage, tirage2)
tirageOn a cette fois une moyenne de 13.312.
On répète l’opération encore 498 fois. On a donc tiré 500 échantillons constitués de 500 personnes. Pour éviter de taper ou exécuter le script 498 fois, on la produit dans une boucle : soit un i allant de 2 à 500, pour chacune des valeurs, on reproduit un tirage aléatoire nouveau et on ajoute la ligne i au tableau des tirages avec le numéro de l’échantillon, la moyenne calculée, l’écart-type calculé et le nombre d’individus dans l’échantillon.
for (i in 2:500) # "éléments" sur lesquels va tourner la boucle
{ # ouverture de la boucle
dt$sample <- sample(1:10000,10000,replace=F) # ajout d'un nombre aléatoire entre 1 et 10000 sans remise
tirage[i,] <- c(i, # Ajout de la ligne i avec le nombre i
mean(dt[dt$sample <501,"choco"]), # la moyenne du tirage de 500
sd = sd(dt[dt$sample <501,"choco"]), # l'écart-type du tirage de 500
n=nrow(dt[dt$sample <501,])) # le nombre de lignes dans le tirage
}
head(tirage)Histogramme de l’ensemble des moyennes de ces 500 échantillons :
ggplot(tirage)+
aes(x=moyenne, stat="count") +
geom_histogram( bins = 20)Sur la base d’un tirage aléatoire simple, en fonction de l’échantillon choisi, la moyenne de la quantité de chocolat mangé l’année dernière varie (ici, assez peu). Cette variabilité est due au hasard de l’échantillon ; on observe une marge d’erreur.
Les conclusions sur la population comportent toujours une marge d’incertitude. Pour être sûr·es à 100% qu’un résultat soit vrai, il faut un échantillon exhaustif ; ie interroger toute la population des 10 000 de départ. Ce n’est que rarement possible. On cherche à probabiliser l’incertain.
On considère, en s’appuyant sur la loi des grands nombres, qu’on s’approche des caractéristiques de la population en augmentant la taille de l’échantillon (tout en ayant l’hypothèse nécessaire de l’aléatoire de l’échantillon). On calcule une moyenne, une proportion et on peut en déduire la probabilité qu’un intervalle contiennent la moyenne ou la proportion pour l’ensemble de la population :
Cette probabilité correspond au seuil de confiance qu’on décide de se donner. L’intervalle contenant – probablement – la moyenne pour la population est l’intervalle de confiance. La valeur calculée sur l’échantillon est finalement une estimation. Note : Il n’est pas sûr que la vraie valeur se trouve dans l’intervalle, ce n’est pas une certitude absolue.
Cette estimation sera en fait un intervalle : On dira qu’« On estime qu’il y a une probabilité de 95% que la moyenne de chocolat consommé chaque année par la population française se situe dans l’intervalle [12,7 ; 13,6].»
L’intervalle de confiance de la moyenne calculée peut être estimé grâce à l’écart-type (dans l’échantillon, la dispersion de la quantité de chocolat consommée l’an dernier est-elle élevée ou non ?). L’estimation à 95% de la moyenne sera l’intervalle :
[ x̅ - (1,96 σ / √n) ; x̅ + (1,96 σ / √n ) ]
avec x̅ : Moyenne ; σ : écart-type ; n : taille de l’échantillon.
Pour le 1er échantillon tiré : Si notre écart-type est de 5.16 l’estimation de la moyenne sera : [ 13.15 – (1.96 x 5.16 / √500) ; 13.15 + (1.96 x 5.16 / √500) ] = [13.15 - ( 1.96 x 5.16 / 22.36) ; 13.15 + ( 1.96 x 5.16 / 22.36) ] = [ 13.15 – 0.45 ; 13.15 + 0.45] = [ 12.7 ; 13.6 ]
On peut dire qu’on a 95 % de chance que la quantité de chocolat consommé soit situé entre :
tirage$moyenne[i] - (1.96 * tirage$sd[i] / (tirage$n[i]^.5) )## [1] 12.34567
et
tirage$moyenne[i] + (1.96 * tirage$sd[i] / (tirage$n[i]^.5) )## [1] 13.20233
La marge d’erreur étant de :
(1.96 * tirage$sd[i] / (tirage$n[i]^.5) )## [1] 0.4283278
et la moyenne calculée de :
mean(tirage$moyenne[i])## [1] 12.774
i<-501
dt$sample <- sample(1:10000,10000,replace=F)
tirage[i,] <- c(i,
mean(dt[dt$sample <51,"choco"]),
sd = sd(dt[dt$sample <51,"choco"]),
n=nrow(dt[dt$sample <51,]))
tirage$moyenne[i] - (1.96 * tirage$sd[i] / (tirage$n[i]^.5) )## [1] 11.6418
tirage$moyenne[i] + (1.96 * tirage$sd[i] / (tirage$n[i]^.5) )## [1] 14.2382
(1.96 * tirage$sd[i] / (tirage$n[i]^.5) )## [1] 1.298198
mean(dt$choco)## [1] 13.1801
L’intervalle de confiance (et donc la précision de l’estimation) est influencé par :
Le niveau de confiance -> au choix, arbitraire, se reporter ensuite à la table de la loi correspondante (ici, loi normale)
La taille de l’échantillon (racine de n) -> suivant les moyens de l’enquête ; plus l’échantillon est grand, plus l’intervalle de confiance est petit
La dispersion des valeurs (écart-type) -> indépendant du chercheur, nécessité de l’observer. Plus il est faible, plus l’intervalle est petit.
Différence entre des échantillons aléatoires constitués de 50 et de 500 individus :
for (i in 502:1000){
dt$sample <- sample(1:10000,10000,replace=F)
tirage[i,] <- c(i,
mean(dt[dt$sample <51,"choco"]),
sd = sd(dt[dt$sample <51,"choco"]),
n=nrow(dt[dt$sample <51,]))
}
ggplot(tirage)+aes(x=moyenne,stat="count") +
geom_histogram() +
stat_bin(bins=30)+
facet_grid(~n)Donc, sur un échantillon donné, les calculs qui seront produits seront assortis d’une marge d’erreur. Il s’agit d’une estimation car l’échantillon reste un échantillon parmi d’autres possibles.
En sciences sociales, on traite surtout des données nominales. Et donc on recherche des proportions (quelle proportion des enquêté·es aiment le chocolat). Prenons donc un deuxième exemple : L’institut de sondage OPIF s’intéresse au plat du dimanche préféré des habitant·es d’un pays. Après un 1er tour élargi, deux plats-candidats restent en lice pour le 2nd : le Dahl de lentilles et le Boeuf bourguignon.
dt2 <- data.frame(repas_pref= c(rep("Dahl de lentille",5100), rep("Boeuf bourguignon", 4900)))La proportion réelle de vote pour le Dahl de lentilles est de 51%. On veut l’estimer en prenant un échantillon de cet ensemble. En reproduisant la démarche du 1er exemple, on tire 1.000 échantillons différents pour voir les variations : le Dahl de lentilles va-t-il l’emporter ?
tirage <- NULL
for (i in 1:1000){
dt2$sample <- sample(1:10000,10000,replace=F)
tirage <- rbind(tirage,
data.frame(tirage=i,
prop=round(nrow(dt2[dt2$sample<501 & dt2$repas_pref=="Dahl de lentille",])/
nrow(dt2[dt2$sample<501,]),3)))
}Graphique des estimations pour ces différents tirages :
ggplot(tirage)+aes(x=prop)+ geom_histogram()+geom_vline(xintercept=.5, col="red")A priori, plus de tirages semblent indiquer que le Dahl l’emporterait mais pas tous. En estimant les résultats du Dahl de lentilles avec un seuil de confiance de 95%, on obtient l’intervalle de confiance suivant :
L’intervalle de confiance d’une proportion s’appuie sur la proportion calcullée et la taille de l’échantillon. L’estimation à 95% de la proportion sera l’intervalle :
[ P - 1.96 * √(P (1 - P) / n) ; P + 1.96 * √(P (1 - P) / n)]
avec P : Probabilité calculée ; n : taille de l’échantillon.
Pour le 1er échantillon tiré : Si notre proportion calculée est de 0.52 l’estimation de la proportion sera :
p=.52
ic= 1.96*(((p*(1-p))/500)^.5)
print(paste0("[",round(p-ic,3)," - ",round(p+ic,3),"]"))## [1] "[0.476 - 0.564]"
Si on souhaite estimer cette proportion avec une probabilité de 90% et qu’on interroge 2.000 personnes :
p=.52
ic= 1.64*(((p*(1-p))/2000)^.5)
print(paste0("[",round(p-ic,3)," - ",round(p+ic,3),"]"))## [1] "[0.502 - 0.538]"
On pourrait alors dire que le Dahl a toutes les chances de l’emporter (avec une probabilité de 90%).
Les tests statistiques permettent d’établir si les résultats obtenus sur l’échantillon sont dus au hasard ou s’ils se rapprochent du résultat qu’on aurait observé sur la population. Ils sonts mobilisés pour confirmer la probabilité (ou plausibilité) d’une corrélation ou d’un lien entre deux variables. La plupart du temps, ils permettent de calculer le risque de se tromper en établissant une hypothèse et ont pour objectif de connaître le niveau de signification (fiabilité) d’un résultat.
On produit une hypothèse appelée l’hypothèse nulle (H0) qui va supposer que les résultats obtenus sont dus au hasard et donc les variables observées sont indépendantes. Le test vise à savoir si on peut rejeter H0 pour affirmer que H1 (hypothèse inverse à H0) est probablement vraie et que les deux variables sont liées.
Pour réaliser et valide ou non ce type d’hypothèses, on va :
A partir d’un tableau d’effectifs croisés, si les proportions semblent montrer que les 2 variables sont liées (jouent dans le même sens ou inversement), le test d’indépendance du khi-2 permet de tester s’il est probable ou non que cette observation soit due au hasard.
Le khi-deux, chi-carré ou chi-square, mesure les écarts entre des effectifs observés et les effectifs théoriques (distribution aléatoire) à partir des marges du tableau croisé (appelé tableau de contingence). Le test du khi2 répond à la question : « Quelle est la probabilité pour que les deux variables qualitatives soient indépendantes l’une de l’autre ? Soit pour que leur distribution se fasse au hasard et qu’il n’y ait aucune relation statistique entre elles ». Il permet d’établir à quel seuil de probabilité on peut rejeter l’hypothèse d’indépendance entre des effectifs observés et les effectifs théoriques. Il nous donne une information sur la corrélation entre deux variables qualitatives.
Un test de khi deux s’applique uniquement sur des tableaux de contingence :
addmargins(table(coco.activite$pcs, coco.activite$sitpro.conf, useNA="ifany"))##
## En congé Au chômage En emploi <NA> Sum
## Agriculteurs 1 2 4 1 8
## Artisans, comm., chef d'E 3 4 12 0 19
## Cadres et pr.intell.sup. 5 14 59 1 79
## Employés 17 20 72 2 111
## Prof. intermédiaires 21 26 55 3 105
## Ouvriers 10 20 38 1 69
## <NA> 21 51 124 3 199
## Sum 78 137 364 11 590
rprop(table(coco.activite$pcs, coco.activite$sitpro.conf, useNA="ifany"))##
## En congé Au chômage En emploi <NA> Total
## Agriculteurs 12.5 25.0 50.0 12.5 100.0
## Artisans, comm., chef d'E 15.8 21.1 63.2 0.0 100.0
## Cadres et pr.intell.sup. 6.3 17.7 74.7 1.3 100.0
## Employés 15.3 18.0 64.9 1.8 100.0
## Prof. intermédiaires 20.0 24.8 52.4 2.9 100.0
## Ouvriers 14.5 29.0 55.1 1.4 100.0
## <NA> 10.6 25.6 62.3 1.5 100.0
## Ensemble 13.2 23.2 61.7 1.9 100.0
Les effectifs sont faibles pour les agriculteur·rices et les artisans. On choisit de les exclure, ainsi que les enquêté·es n’ayant pas donné d’élément permettant de connaître leur situation professionnelle pendant le confinement.
coco.activite.khi <- coco.activite %>%
filter(! pcs %in% c("Agriculteurs",
"Artisans, comm., chef d'E")) %>%
filter(is.na(sitpro.conf)==F) %>%
filter(is.na(pcs)==F) %>%
mutate(pcs = fct_drop(pcs))
addmargins(table(coco.activite.khi$pcs, coco.activite.khi$sitpro.conf, useNA="ifany"))##
## En congé Au chômage En emploi Sum
## Cadres et pr.intell.sup. 5 14 59 78
## Employés 17 20 72 109
## Prof. intermédiaires 21 26 55 102
## Ouvriers 10 20 38 68
## Sum 53 80 224 357
rprop(table(coco.activite.khi$pcs, coco.activite.khi$sitpro.conf, useNA="ifany"))##
## En congé Au chômage En emploi Total
## Cadres et pr.intell.sup. 6.4 17.9 75.6 100.0
## Employés 15.6 18.3 66.1 100.0
## Prof. intermédiaires 20.6 25.5 53.9 100.0
## Ouvriers 14.7 29.4 55.9 100.0
## Ensemble 14.8 22.4 62.7 100.0
coco.activite.khi %>%
ggplot()+
aes(fill=sitpro.conf, x = pcs)+
geom_bar(position="fill")+
labs(title = "Situation pro. pendant le confinement des personnes en activité avant selon la PCS",
x="PCS",fill = "Situation professionnelle",y = "Proportion")A première vue, les cadres et professions intellectuelles supérieures sont plus nombreux·ses à être resté·es en emploi, quand les professions intermédiaires ont davantage été mises en congés et les ouvriers et ouvrières plus souvent été au chômage pendant cette période. Il s’agit maintenant de vérifier ce résultat avec un test du khi-deux :
chisq.test(table(coco.activite.khi$pcs, coco.activite.khi$sitpro.conf))##
## Pearson's Chi-squared test
##
## data: table(coco.activite.khi$pcs, coco.activite.khi$sitpro.conf)
## X-squared = 13.495, df = 6, p-value = 0.03582
Cette instruction donne le résultat pour l’ensemble du tableau :
Par convention, les seuils suivants sont retenus :
Pour l’exemple, l’association observée est significative. On peut affirmer qu’il est probable que la pcs soit liée à la situation professionnelle durant le confinement des personnes qui étaient auparavant en activité.
…Une dernière chose :
# Produire le survey design :
coco.activite.khi.w <- svydesign(data=coco.activite.khi,
ids=~1, weights = coco.activite.khi$POIDS)
# Tri croisé :
svytable(~ pcs + sitpro.conf, coco.activite.khi.w) %>%
rprop(2)## sitpro.conf
## pcs En congé Au chômage En emploi Total
## Cadres et pr.intell.sup. 11.40 19.56 69.04 100.00
## Employés 14.29 17.29 68.42 100.00
## Prof. intermédiaires 23.93 28.94 47.12 100.00
## Ouvriers 28.66 26.71 44.64 100.00
## Ensemble 20.23 23.74 56.03 100.00
svytable(~ pcs + sitpro.conf, coco.activite.khi.w) %>%
addmargins() %>% round(0)## sitpro.conf
## pcs En congé Au chômage En emploi Sum
## Cadres et pr.intell.sup. 7 12 43 62
## Employés 13 15 60 88
## Prof. intermédiaires 28 33 54 115
## Ouvriers 22 21 34 77
## Sum 69 81 192 343
# Test du chi-2 pondéré :
svychisq(~ pcs + sitpro.conf, coco.activite.khi.w, statistic="Chisq")##
## Pearson's X^2: Rao & Scott adjustment
##
## data: svychisq(~pcs + sitpro.conf, coco.activite.khi.w, statistic = "Chisq")
## X-squared = 19.516, df = 6, p-value = 0.2267
Cette fois-ci, le résultat tend à démontrer qu’il y a une forte probabilité pour que le lien visible entre les 2 variables soit du au hasard. Du fait des filtres que nous avons du faire, l’échantillon final sur lequel porte le test statistique est très faible.
Attention :
Plus de détails : Voir Barnier Julien, « Tout ce que vous n’avez jamais voulu savoir sur le χ² sans jamais avoir eu envie de le demander », 2016
Si le test du khi-2 est l’un des plus utilisés en sciences sociales (car il porte sur les fréquences), il existe une multitude de tests statistiques selon les configurations des données. Chacun porte ses propres conditions et hypothèses.
Voir aussi : https://larmarange.github.io/analyse-R/comparaisons-moyennes-et-proportions.html
Pour les estimations : https://larmarange.github.io/analyse-R/intervalles-de-confiance.html
Le test du Khi-Deux a une autre limite : il permet de tester l’indépendance entre deux variables nominales, mais comment fait-on pour mettre à l’épreuve la relation entre une variable qualitative et une variable quantitative (une variable nominale et une variable numérique) ? Parmi les tests possibles, le moins mal connu et le moins inutilisé en sciences sociales est l’ANOVA pour analyse de variance qui compare les variances entre différents groupes.
On compare la variance dans les groupes et entre ces groupes (intergroupe et intragroupe), les différents groupes correspondant aux modalités d’une variable catégorielle :
Si la variance intragroupe est supérieure à la variance intergroupe, on peut en conclure que le groupe d’appartenance ne modifie pas les valeurs de la variable numérique prise en compte : il n’y a pas de lien entre les variables ;
Si la variance intergroupe est supérieure à la variance intragroupe, alors c’est le signe d’une relation entre les deux variables ;
Sur le même modèle que le Khi-2, l’ANOVA compare les variances observées avec les variances théoriques en cas d’indépendance totale entre les deux variables ; l’hypothèse nulle est donc, comme pour le Khi-Deux, l’indépendance ;
Aux limites présentées pour le Khi-Deux s’en ajoute une autre : la variable numérique doit réunir certaines propriétés, notamment la normalité (avoir une distribution suivant une loi normale) ;
Dans R :
anova <- aov(variableNumerique ~ variableNominale)
Où le ~ s’écrit AltGR+2 sous un PC, Alt enfoncé + N sous Mac) ;
Puis summary(anova)pour obtenir le
résultat.
Les fichiers et la synthaxe Rmarkdown permettent de
générer des rapports automatiques mêlant du texte, du
code et les résultats sur des jeux de
données. Ce fichier contient :
Si vous êtes déjà dans un projet R (.Rproj), vous pouvez créer un fichier .Rmd :
Fichier > Nouveau fichier > Fichier R Markdown ; puis ok dans la fenêtre. Un fichier avec du texte d’exemple par défaut apparaît. Pour exécuter l’ensemble et produire le document : bouton “knit” ou ctrl + alt + K et nommer le fichier en sortie (uniquement la 1ère fois).
Exemple : Les présentations supports du cours ont été générées avec
rmarkdown
L’en-tête est obligatoire et définit les principales caractéristiques du document produit. L’en-tête par défaut est minimale. Elle est donc définie en tout début de document et est délimitée par deux lignes “—”.
Attention, cet espace respecte l’indentation.
output: html_document()ou
output: pdf_document()ou
output: word_document()Voir : https://bookdown.org/yihui/rmarkdown/output-formats.html
title: "Titre du document"author: "Nom de l'auteur"date: 'Dernière modification: `r format(Sys.Date(), format="%d/%m/%Y")`'On peut voir que des instructions R apparaissent. Elles sont balisées comme ceci :
`r code R `subtitle: "Sous-titre"abstract: "Ajout d'un résumé\nRetour à la ligne"Le format généré permet d’intégrer une série d’options dont :
number_sections: yes # pour ajouter des numéros aux titres
toc: yes # ajout d'une table des matières
toc_depth : '3' # Précision de la table des matières
folding_code : hide # Cacher le code par défaut
toc_float: true # Laisser la table des matières constamment visibleEt une série d’autres options précisées dans la cheat sheet Rmarkdown : https://raw.githubusercontent.com/rstudio/cheatsheets/master/rmarkdown-2.0.pdf
---
title: "R - Avancé - M2 ENSP"
author: "Cécile Rodrigues - CNRS / CERAPS"
date: 'Dernière modification : `r format(Sys.Date(), format="%d/%m/%Y")`'
output:
html_document:
number_sections: yes
toc: yes
toc_depth: '3'
toc_float: true
---Rappel : L’indentation compte !
Ici, on a ajouté des numéros aux titres de sections, une table des matières allant jusqu’au 3ème niveau et présenté sur le côté (au lieu d’être en début de document).
Les titres sont précédés par un # pour les titres de niveau 1 ; par ## pour les titres de niveau 2, etc.
Le texte brut écrit dans le fichier .rmd apparaîtra comme le texte par défaut d’un traitement de texte. Le saut de paragraphe se fait en sautant deux lignes.
Les tirets seront matérialisés par des puces :
- Comme
- ici
- En ajoutant des sous-catégories
1. Ou bien
2. avec des
3. numérosLa mise en forme italique est indiquée en entourant le texte à metttre en italique par * : *texte en italique* ; et celui en gras par ** (** texte en gras **)
Du code peut être inséré directement dans le texte :
`r 1+1 `Le code et les résultats R sont directement intégrés dans le texte mais balisé par
Les balises apparaissent en cliquant sur insert > R / python / etc ou en tapant Ctrl+Alt+I.
Exemple :
Sys.Date()Donne :
## [1] "2022-11-23"
Plusieurs options peuvent être ajoutées dans l’en-tête du chunk :
{ r nom du chunk, eval = F, echo = F, warning=F }
Détail des options :
Pour définir ces options systématiquement, ajouter en début de script .Rmd :
knitr::opts_chunk$set(warning=FALSE, eval=FALSE, message=FALSE, error=FALSE)Les tableaux et graphiques sont générés et mis en forme dans les chunks. Plusieurs packages permettent de produire des tableaux esthétiques, parfois complexes, et contenant des mises en forme particulières (couleur de fond, texte en gras, ajout d’une barre en fonction d’une valeur, fusionner des cellules, etc.). Dans le détail, ces fonctions produisent des tableaux en code html ou latex (suivant la fonction et le format de sortie) :
un premier package utile pour faire des tableaux est {knitr} et les
fonctions kable. Dans sa forme la plus basique, en
reprenant l’exemple du tableau des “vies communes” :
vitavec <- freqm(coco, "q24", calc="%", exclude="6666", poids="POIDS") %>%
mutate(name= fct_recode(name,
"Conjoint·e"="q24_1",
"Enfants/petits-enfants"="q24_2",
"Parents/Grands-parents"="q24_3",
"Famille (autre)"="q24_4",
"Ami·es ou autres"="q24_5",
"Chien(s)"="q24_6",
"Autre(s) animal(aux)"="q24_7")) %>%
rename("Vit avec"="name", "Oui"="X1") %>%
select(-X2, -Total) %>% arrange(desc(Oui))library(knitr)
vitavec %>%
kable("html", escape = F, caption = "Source : Coco1 - ELIPSS, CDSP, avril 2020")| Vit avec | Oui |
|---|---|
| Conjoint·e | 81.6 |
| Enfants/petits-enfants | 48.5 |
| Autre(s) animal(aux) | 23.5 |
| Chien(s) | 14.2 |
| Parents/Grands-parents | 7.1 |
| Famille (autre) | 5.6 |
| Ami·es ou autres | 3.7 |
En un peu plus stylisé (en version html_document) :
library(kableExtra)
vitavec %>%
kable("html", escape = F, caption = "Source : Coco1 - ELIPSS, CDSP, avril 2020") %>%
kable_styling("striped", full_width = F)| Vit avec | Oui |
|---|---|
| Conjoint·e | 81.6 |
| Enfants/petits-enfants | 48.5 |
| Autre(s) animal(aux) | 23.5 |
| Chien(s) | 14.2 |
| Parents/Grands-parents | 7.1 |
| Famille (autre) | 5.6 |
| Ami·es ou autres | 3.7 |
-> Améliorer les tableaux avec {kableExtra} : https://cran.r-project.org/web/packages/kableExtra/vignettes/awesome_table_in_html.html
Les graphiques de type ggplot deviennent des images mais vous pouvez agir sur leur taille. Ajouter les dimensions du graphique dans le chunk :
{r, fig.width= 9, fig.height=5}
vitavec %>% arrange(desc(Oui)) %>% mutate(`Vit avec` = fct_inorder(`Vit avec`)) %>%
ggplot() + geom_histogram(aes(x = `Vit avec`, y = Oui), stat="identity", col="black",fill="white") +
theme_bw()+ theme(axis.text.x = element_text(angle=45, hjust=1))+
labs(y="Pourcentages (%)",
title="Avec qui vivent les enquêté·es", caption="Source : Coco1 - ELIPSS, CDSP, avril 2020")2 autres packages intéressants :
Les graphiques de type ggplot deviennent des images mais vous pouvez agir sur leur taille. Ajouter les dimensions du graphique dans le chunk :
{r, fig.width= 9, fig.height=5}
vitavec %>% arrange(desc(Oui)) %>% mutate(`Vit avec` = fct_inorder(`Vit avec`)) %>%
ggplot() + geom_histogram(aes(x = `Vit avec`, y = Oui), stat="identity", col="black",fill="white") +
theme_bw()+ theme(axis.text.x = element_text(angle=45, hjust=1))+
labs(y="Pourcentages (%)",
title="Avec qui vivent les enquêté·es", caption="Source : Coco1 - ELIPSS, CDSP, avril 2020")Enfin, les images peuvent être ajoutées de deux façons :
En ajoutant dans l’en-tête out.width = ‘400px’ pour préciser la taille :
knitr::include_graphics("img/tidydata_7.jpg")Cet été, la société qui produit RStudio a fait deux annonces. Elle s’appelle désormais posit mais le nom du logiciel reste rstudio et elle a lancé un nouveau projet : les fichiers quarto. Quarto est pensé comme un nouvel outil de publication scientifique et technique se rapprochant de rmarkown et des notebook (jupyter, par exemple). L’idée est de produire des documents hybrides, contenant du texte, de la rédaction, des commentaires et analyses mais aussi des figures, tables et résultats directement issus de codes intégrés dans le document. Quarto est présenté comme la “nouvelle génération” de rmarkdown et pousse l’outil sur d’autres langages (python, julia, javascript, etc.) ; ce qui était possible en rmarkdown sans que ça soit au centre de l’outil.
Voir : https://posit.co/blog/announcing-quarto-a-new-scientific-and-technical-publishing-system/
Un exemple : https://neocarto.github.io/docs/notebooks/bertinquarto/
–
Pour aller plus loin dans l’interactivité, vous pourrez utiliser le langage shiny qui permet de créer des applications web interactives dans l’éco-système de R. Cf l’intervention que vous aurez sur ce langage.
ctrl+alt+I : créer un chunk
ctrl+shift+K : “knit”, produire le document
Comme on est passé·es un peu rapidement sur les exercices, prenons le temps de les faire maintenant :
Dans R, on crée des objets et on leur applique des fonctions pour produire des résultats et des analyses. Les fonctions prennent (1) un ou plusieurs arguments en entrée (ou paramètres ou input), éventuellement certaines options, (2) effectuent des opérations ou des actions et (3) produisent et renvoient des résultats ou sorties (output).
R est quand même un langage de programmation qui permet de créer ses propres fonctions si on ne les trouve pas dans le langage de base ou via un package.
L’écriture, dans sa version la plus basique, d’une fonction sera :
Ajoute2 <- function(x){
res <- x + 2
return(res)
}
Ajoute2(4)objet <-.Ajoute2 est apparu dans la partie “Functions” de
l’environnement.
Attention : Si le nom de la fonction existait déjà, l’ancienne fonction est écrasée
summary <- function(x){
res <- x+2
return(res)}summary est maintenant remplacée par le contenu de
Ajoute2.
Mais on peut toujours supprimer la fonction summary
qu’on vient de créer localement :
rm(summary)propTable <- function(v){
tri <- table(v)
total <- length(v)
tri <- round(tri/total*100,2)
return(tri)
}
propTable(coco$q01)## v
## 1 2
## 19.89 80.11
Ca peut être une fonction qui intègre la création d’un graphique :
MyBarplot <- function(var){
tri <- table(var)
barplot(tri, col ="skyblue", border=NA)
}
MyBarplot(coco$q01)Les fonctions sont très utiles quand on se retrouve à répéter des actions. Par exemple, on recode une dizaine de variables sur le même modèle : 1 = Oui et 2 = Non :
Au lieu d’écrire…
library(questionr)
library(tidyverse)
freq(coco$q24_1)coco$q24_1_r <- fct_collapse(factor(coco$q24_1),
"Oui"="1", "Non"="2", "NR"="6666")
freq(coco$q24_1_r)… On crée la fonction qui fait ce recodage et on l’applique :
fct_ON <- function(x){
recod <- fct_collapse(factor(x),
"Oui"="1", "Non"="2", "NR"="6666")
return(recod)
}
coco$q24_2_r <- fct_ON(coco$q24_2)
freq(coco$q24_2_r)Les intérêts sont : - On évite une série de copier/coller qui alourdissent le code ; - Plus simple pour corriger d’éventuelles erreurs : on corrige la fonction et pas les 7 copies - Réutilisable d’un script à l’autre, vous pouvez faire des scripts avec vos fonctions habituelles
Comme pour les fonctions en général, les arguments peuvent être
explicités ou être implicites selon leur position :
fct_ON(x = coco$q24_2) ou fct_ON(coco$s24_2)
car il n’y a pour l’instant qu’un seul argument.
En général, on compte sur l’implicite pour les données et on
explicite les options; Comme dans mean(x, na.rm = T)
Un autre intérêt sera de tester un morceau de code en changeant un ou quelques paramètres :
graphede <- function(x){
des <- data.frame(jets=sample(1:6, x, rep=T))
ggplot(des)+aes(x=jets)+geom_bar(fill="white",col="black")+theme_bw()+scale_x_continuous(breaks = 1:6)
}
graphede(500)Il est possible d’intégrer des valeurs par défaut comme paramètres :
Nom <- function(x, dec=1, useNA="ifany"){}–
Sans valeur par défaut, un argument est obligatoire mais facultatif si une valeur par défaut apparaît.
En ajoutant …, on indique qu’on intègre les arguments par défaut des fonctions appelées dans celle qu’on crée :
Nom <- function(x, ...){}Par défaut, R renvoie le dernier objet présent dans la fonction. Mais
en utilisant la fonction return, on peut davantage
maîtriser ce qui est renvoyé. C’est par exemple utile quand on veut
renvoyer un élément différent sur la base de conditions.
Comportement pas défaut :
fun <- function(x){
x+2
x+5
}
fun(2)–
Pour forcer le retour de la 1ère ligne :
fun <- function(x){
return(x+2)
x+5}
fun(2)Créez une fonction qui renvoie le carré et le cube d’un nombre : ^
pour élever un chiffre à une puissance données : 2^2.
–
carrecube <- function(x) {
carre = x^2
cube = x^3
return(list(x = x,
carre = carre,
cube=cube))
}On renvoie une liste qui contient le nombre mis en entrée, son carré
et son cube. Si on fait carrecube("Coucou") ?
–
On obtient une erreur qu’on peut éviter avec return car
une fois qu’une fonction return est actionnée, le reste de
la fonction n’est pas exécuté.
carrecube <- function(x) {
if(!is.numeric(x)) return("Eh non. Il faut mettre un chiffre !")
carre = x^2
cube = x^3
return(list(x = x,
carre = carre,
cube=cube))}
carrecube("coucou")## [1] "Eh non. Il faut mettre un chiffre !"
Note : Si on met “print” au lieu de “return”, le code continue à s’exécuter et on rencontre une erreur.
On n’est pas obligé·es de définir les arguments dans la fonction (local), ils peuvent être dans l’environnement global (la session rstudio). Mais si le même objet est dans la fonction, il sera priorisé.
Il est aussi possible de modifier un objet dans l’environnement
global à partir d’une fonction en utilissant l’opérateur
<<- :
x <- 4
changex <- function(a){x <<- a}
x## [1] 4
changex(5)
x## [1] 5
Un nombre est premier s’il n’est divisible que par lui-même et par 1, c’est-à-dire si le reste de la division est égal à 0 uniquement quand on le divise par lui-même ou 1.
4%%2premier <- function(x){
# Pas tout seul-----
if(length(x)>1) return(print("x doit être un nombre seul"))
# Pas un nombre --------
if(! class(x) %in% c("integer", "numeric")) return(print("x doit être un nombre"))
# Pas entier ---------
if(grepl("\\.",x)) return(print("x doit être un nombre entier"))
# le nombre est-il premier ? -------
for(i in (2:(x^.5))){
if(x %% i == 0) return(paste0(x, " n'est pas premier"))}
return(paste0(x, " est premier"))
}
premier(53)## [1] "53 est premier"
listepremier <- function(x,y){
liste <- NULL
for (i in (x:y)){
for(j in (2:(i^.5))){
if(i %% j == 0){liste <- c(liste, i)}
}}
premiers <- setdiff(x:y, liste)
return(paste0(paste0("Les nombres premiers entre ",x," et ",y, " sont : "),
paste(premiers, collapse = ", ")))
}
listepremier(10,20)## [1] "Les nombres premiers entre 10 et 20 sont : 11, 13, 17, 19"
Il s’agira d’un dossier en groupe de 2 à 4 étudiant·es au choix à partir de données fournies et à rechercher autour d’une même thématique. Le dossier rendu sera :
-> Vous m’enverrez le document produit et le script .Rmd au plus tard le 03/01/23 à mon adresse mail cecile.rodrigues@univ-lille.fr (vous aurez un accusé de réception sous 24h).
Source : allison
Horst @allison_horst
#rstats
#rstatsFR
@icymi_r
@hadleywickham
@thinkR_fr
@Rbloggers
@tidyversetweets